fn_rolling_stocks_par <- function(df_input, controls_levels, controls_logs, window_T, unlearnable, quarter_lag, group_type="permno"){

  group_list <- df_input %>% group_by(across(all_of(group_type))) %>% count()
  n_groups    <- group_list %>% nrow() %>% as.numeric()
  
  n_outputs <- 17
  
  df <- df_input
  
  controls_logs_lagged <- c()
  controls_levels_lagged <- c()
  
  if (length(controls_logs) != 0 ) {
    controls_logs_lagged <- vector(length = length(controls_logs))
    for (i in 1:length(controls_logs)) {
      varname <- paste("delta_log_", controls_logs[i], sep = "")
      controls_logs_lagged[i] = varname
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname := log(get(controls_logs[i])) - log(dplyr::lag(get(controls_logs[i]), n = quarter_lag, default = NA))) %>%
        ungroup()
    }
  }
  
  if (length(controls_levels) != 0 ) {
    controls_levels_lagged <- vector(length = length(controls_levels))
    for (i in 1:length(controls_levels)) {
      varname <- paste("delta_", controls_levels[i], sep = "")
      controls_levels_lagged[i] = varname
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname := get(controls_levels[i]) - dplyr::lag(get(controls_levels[i]), n = quarter_lag, default = NA)) %>%
        ungroup()
    }
  }
  
  dep_var <- "delta_log_price"
  
  if (length(controls_levels_lagged) != 0 | length(controls_logs_lagged) != 0) {
    controls_lagged = append(controls_logs_lagged, controls_levels_lagged)
  } else {
    controls_lagged = c()
  }
  
  if (unlearnable) {
    indep_var_long  <- c("delta_log_payoff_growth", "forecast_growth", "future_forecast_growth")
    indep_var_short <- c("delta_log_payoff_growth", "forecast_growth")
  } else {
    indep_var_long  <- c("delta_log_payoff_growth", "delta_log_payoff_growth_future")
    indep_var_short <- c("delta_log_payoff_growth")
  }
  
  indep_var_long  <- append(indep_var_long, controls_lagged)
  indep_var_short <- append(indep_var_short, controls_lagged)
  
  pb <- txtProgressBar(min=1, max=n_groups, style=3)
  progress <- function(n) setTxtProgressBar(pb, n)
  opts <- list(progress=progress)
  
  results <- foreach(i = 1:n_groups, .combine = "rbind", .export = c('fn_identifying_stocks'), .packages = 'tidyverse', .options.snow=opts) %dopar% {
    group       <- group_list[i,1] %>% as.numeric()
    
    df_group <- df %>% filter(.[[group_type]] == group)
    
    T_periods <- df_group %>% count() %>% pull() %>% as.numeric() # Total observations
    T_rolling <- T_periods - window_T + 1        # Number of iterations
    
    output_i <- matrix(ncol = n_outputs, nrow = T_rolling) # Empty data frame to store results
    
    for (j in 1:T_rolling) {
      df_window <- df_group %>% slice(j:(j + window_T - 1)) ## Selecting relevant subsample
      
      year_start  <- year(min(df_window$date_fred))
      month_start <- month(min(df_window$date_fred))
      year_end    <- year(max(df_window$date_fred))
      month_end   <- month(max(df_window$date_fred))
      
      df_window <- df_window %>% filter(if_all(append(indep_var_long, dep_var), is.finite))
      
      recovered_list <- fn_identifying_stocks(df_window, group, dep_var, indep_var_long, indep_var_short)
      
      output_i[j, 1:4] <- c(year_start, month_start, year_end, month_end)
      output_i[j, 5:n_outputs] <- unlist(recovered_list, use.names = FALSE)
      
    }
    output_i
  }
  
  close(pb)
  
  output <- data.frame(results)
  colnames(output) <- c("year_start", "month_start", "year_end", "month_end", 
                        group_type, "n_obs", "n_factors", "R2_long", "R2_short",
                        "beta_0", "beta_1", "zeta_0", "var_eps", "out_long", "out_short", 
                        "tau_pi", "tau_pi_R")
  
  return(output)
  
}

fn_rolling_stocks_new_par <- function(df_input, controls_levels, controls_logs, window_T, unlearnable, quarter_lag, group_type="permno"){
  
  group_list <- df_input %>% group_by(across(all_of(group_type))) %>% count()
  n_groups    <- group_list %>% nrow() %>% as.numeric()
  
  n_outputs <- 17
  
  df <- df_input
  
  controls_logs_lagged <- c()
  controls_levels_lagged <- c()
  
  if (length(controls_logs) != 0 ) {
    controls_logs_lagged <- vector(length = length(controls_logs))
    for (i in 1:length(controls_logs)) {
      varname <- paste("delta_log_", controls_logs[i], sep = "")
      controls_logs_lagged[i] = varname
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname := log(get(controls_logs[i])) - log(dplyr::lag(get(controls_logs[i]), n = quarter_lag, default = NA))) %>%
        ungroup()
    }
  }
  
  if (length(controls_levels) != 0 ) {
    controls_levels_lagged <- vector(length = length(controls_levels))
    for (i in 1:length(controls_levels)) {
      varname <- paste("delta_", controls_levels[i], sep = "")
      controls_levels_lagged[i] = varname
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname := get(controls_levels[i]) - dplyr::lag(get(controls_levels[i]), n = quarter_lag, default = NA)) %>%
        ungroup()
    }
  }
  
  dep_var <- "delta_log_price"
  
  if (length(controls_levels_lagged) != 0 | length(controls_logs_lagged) != 0) {
    controls_lagged = append(controls_logs_lagged, controls_levels_lagged)
  } else {
    controls_lagged = c()
  }
  
  if (unlearnable) {
    indep_var_long  <- c("delta_log_payoff_growth", "forecast_growth", "future_forecast_growth")
    indep_var_short <- c("delta_log_payoff_growth", "forecast_growth")
  } else {
    indep_var_long  <- c("delta_log_payoff_growth", "delta_log_payoff_growth_future")
    indep_var_short <- c("delta_log_payoff_growth")
  }
  
  indep_var_long  <- append(indep_var_long, controls_lagged)
  indep_var_short <- append(indep_var_short, controls_lagged)
  
  pb <- txtProgressBar(min=1, max=n_groups, style=3)
  progress <- function(n) setTxtProgressBar(pb, n)
  opts <- list(progress=progress)
  
  results <- foreach(i = 1:n_groups, .combine = "rbind", .export = c('fn_identifying_stocks'), .packages = 'tidyverse', .options.snow=opts) %dopar% {
    group       <- group_list[i,1] %>% as.numeric()
    
    df_group <- df %>% filter(.[[group_type]] == group)
    
    periods <- df_group %>% select(date_fred) %>% unique()
    
    T_periods <- periods %>% count() %>% pull() %>% as.numeric() # Total observations
    T_rolling <- T_periods - window_T + 1        # Number of iterations
    
    output_i <- matrix(ncol = n_outputs, nrow = T_rolling) # Empty data frame to store results
    
    for (j in 1:T_rolling) {
      period_start <- periods$date_fred[j]
      period_end <- periods$date_fred[j + window_T - 1]
      df_window <- df_group %>% filter(between(date_fred, period_start, period_end)) ## Selecting relevant subsample
      
      year_start  <- year(min(df_window$date_fred))
      month_start <- month(min(df_window$date_fred))
      year_end    <- year(max(df_window$date_fred))
      month_end   <- month(max(df_window$date_fred))
      
      df_window <- df_window %>% filter(if_all(append(indep_var_long, dep_var), is.finite))
      
      recovered_list <- fn_identifying_stocks(df_window, group, dep_var, indep_var_long, indep_var_short)
      
      output_i[j, 1:4] <- c(year_start, month_start, year_end, month_end)
      output_i[j, 5:n_outputs] <- unlist(recovered_list, use.names = FALSE)
      
    }
    output_i
  }
  
  close(pb)
  
  output <- data.frame(results)
  colnames(output) <- c("year_start", "month_start", "year_end", "month_end", 
                        group_type, "n_obs", "n_factors", "R2_long", "R2_short",
                        "beta_0", "beta_1", "zeta_0", "var_eps", "out_long", "out_short",
                        "tau_pi", "tau_pi_R")
  
  return(output)
  
}

build_bootstrap_series <- function(df_i, boot_blocks, window_T){
  nrows <- df_i %>% nrow() %>% as.numeric()
  indices <- numeric(nrows)
  w <- 1
  for(i in 1:nrows){
    block_pos <- boot_blocks$block_starts[[i]]
    block_length <- boot_blocks$block_lengths[[i]]
    if(block_pos > nrows){
      block_length <- block_length + nrows - block_pos
      block_pos <- 1
    }
    j <- 0
    for(jj in 1:block_length){
      indices[w] <- block_pos + j
      w <- w + 1
      j <- j + 1
      if(block_pos + j > nrows){
        j <- 1 - block_pos
      }
      if(w > nrows){
        break
      }
    }
    if(w > nrows){
      break
    }
  }
  
  output <- df_i[indices,]
  return(output)
}

fn_rolling_stocks_boot_par <- function(df_input, controls_levels, controls_logs, window_T, unlearnable, quarter_lag, boot_samples=10, block_override=-1, group_type="permno"){
  group_list <- df_input %>% group_by(across(all_of(group_type))) %>% count()
  n_groups    <- group_list %>% nrow() %>% as.numeric()
  
  n_outputs <- 17
  
  df <- df_input
  
  df <- df %>% mutate(quarter = ceiling(month/3))
  
  controls_logs_lagged <- c()
  controls_levels_lagged <- c()
  
  if (length(controls_logs) != 0 ) {
    controls_logs_lagged <- vector(length = length(controls_logs))
    for (i in 1:length(controls_logs)) {
      varname <- paste("delta_log_", controls_logs[i], sep = "")
      controls_logs_lagged[i] = varname
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname := log(get(controls_logs[i])) - log(dplyr::lag(get(controls_logs[i]), n = quarter_lag, default = NA))) %>%
        ungroup()
    }
  }
  
  if (length(controls_levels) != 0 ) {
    controls_levels_lagged <- vector(length = length(controls_levels))
    for (i in 1:length(controls_levels)) {
      varname <- paste("delta_", controls_levels[i], sep = "")
      controls_levels_lagged[i] = varname
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname := get(controls_levels[i]) - dplyr::lag(get(controls_levels[i]), n = quarter_lag, default = NA)) %>%
        ungroup()
    }
  }
  
  dep_var <- "delta_log_price"
  
  if (length(controls_levels_lagged) != 0 | length(controls_logs_lagged) != 0) {
    controls_lagged = append(controls_logs_lagged, controls_levels_lagged)
  } else {
    controls_lagged = c()
  }
  
  if (unlearnable) {
    indep_var_long  <- c("delta_log_payoff_growth", "forecast_growth", "future_forecast_growth")
    indep_var_short <- c("delta_log_payoff_growth", "forecast_growth")
  } else {
    indep_var_long  <- c("delta_log_payoff_growth", "delta_log_payoff_growth_future")
    indep_var_short <- c("delta_log_payoff_growth")
  }
  
  indep_var_long  <- append(indep_var_long, controls_lagged)
  indep_var_short <- append(indep_var_short, controls_lagged)
  
  expected_blocklength <- 1
  if(block_override > 0){
    expected_blocklength <- block_override
  }else{
    expected_blocklength <- round(mean(unlist(df %>% group_by(across(all_of(group_type))) %>% group_map(~pwsd(.x$delta_log_price, correlogram = F)[["BlockLength"]][1]))))
  }
  
  blocks <- crossing(boot_sample=1:boot_samples, block=1:window_T)
  c <- blocks %>% count() %>% as.numeric()
  
  blocks$block_starts <- ceiling(runif(c, min=0, max=window_T))
  blocks$block_lengths <- 1+rgeom(c, 1/(expected_blocklength-1))
  
  pb <- txtProgressBar(min=1, max=n_groups, style=3)
  progress <- function(n) setTxtProgressBar(pb, n)
  opts <- list(progress=progress)
  
  results <- foreach(i = 1:n_groups, .combine = "rbind", .export = c('fn_identifying_stocks', 'build_bootstrap_series'), .packages = 'tidyverse', .options.snow=opts) %dopar% {
    group <- group_list[i,1] %>% as.numeric()
    
    df_group <- df %>% filter(.[[group_type]] == group)
    
    T_periods <- df_group %>% count() %>% pull() %>% as.numeric() # Total observations
    T_rolling <- T_periods - window_T + 1        # Number of iterations
    
    output_i <- matrix(nrow = boot_samples*T_rolling, ncol = n_outputs+1) # Empty data frame to store results
    
    for (j in 1:T_rolling) {
      df_window <- df_group %>% slice(j:(j + window_T - 1)) ## Selecting relevant subsample
      
      year_start  <- year(min(df_window$date_fred))
      month_start <- month(min(df_window$date_fred))
      year_end    <- year(max(df_window$date_fred))
      month_end   <- month(max(df_window$date_fred))
      
      df_window <- df_window %>% filter(if_all(append(indep_var_long, dep_var), is.finite))
      
      for (k in 1:boot_samples){
        curr_blocks <- blocks %>% filter(boot_sample==k)
        
        block_sample <- build_bootstrap_series(df_window, curr_blocks, window_T)
        
        recovered_list <- fn_identifying_stocks(block_sample, group, dep_var, indep_var_long, indep_var_short)
        output_i[boot_samples*(j-1)+k, 1:4] <- c(year_start, month_start, year_end, month_end)
        output_i[boot_samples*(j-1)+k, 5:n_outputs] <- unlist(recovered_list, use.names = FALSE)
        output_i[boot_samples*(j-1)+k, n_outputs+1] <- k
      }
    }
    output_i
  }
  
  close(pb)
  
  output <- data.frame(results)
  colnames(output) <- c("year_start", "month_start", "year_end", "month_end", 
                        group_type, "n_obs", "n_factors", "R2_long", "R2_short",
                        "beta_0", "beta_1", "zeta_0", "var_eps", "out_long", "out_short",
                        "tau_pi", "tau_pi_R", "boot_run")
  
  return(output)
  
}

fn_rolling_stocks_PS_par <- function(df_input, ivol, df_io_share, controls_levels, controls_logs, window_T, unlearnable, quarter_lag, group_type="permno"){
  group_list <- df_input %>% group_by(across(all_of(group_type))) %>% count()
  n_groups    <- group_list %>% nrow() %>% as.numeric()
  
  df <- df_input
  
  controls_logs_lagged <- c()
  controls_levels_lagged <- c()
  controls_logs_inter <- c()
  controls_levels_inter <- c()
  controls_logs_inter2 <- c()
  controls_levels_inter2 <- c()
  
  if (length(controls_logs) != 0 ) {
    controls_logs_lagged <- vector(length = length(controls_logs))
    for (i in 1:length(controls_logs)) {
      varname_lag <- paste("delta_log_", controls_logs[i], sep = "")
      varname_inter <- paste("inter_log_", controls_logs[i], sep = "")
      varname_inter2 <- paste("inter2_log_", controls_logs[i], sep = "")
      controls_logs_lagged[i] = varname_lag
      controls_logs_inter[i] = varname_inter
      controls_logs_inter2[i] = varname_inter2
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname_lag := log(get(controls_logs[i])) - log(dplyr::lag(get(controls_logs[i]), n = quarter_lag, default = NA))) %>%
        mutate(!!varname_inter := get(varname_lag) * delta_log_payoff_growth_future) %>%
        mutate(!!varname_inter2 := get(varname_lag) * delta_log_payoff_growth) %>%
        ungroup()
    }
  }
  
  if (length(controls_levels) != 0 ) {
    controls_levels_lagged <- vector(length = length(controls_levels))
    for (i in 1:length(controls_levels)) {
      varname_lag <- paste("delta_", controls_levels[i], sep = "")
      varname_inter <- paste("inter_", controls_levels[i], sep = "")
      varname_inter2 <- paste("inter2_", controls_levels[i], sep = "")
      controls_levels_lagged[i] = varname_lag
      controls_levels_inter[i] = varname_inter
      controls_levels_inter2[i] = varname_inter2
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname_lag := get(controls_levels[i]) - dplyr::lag(get(controls_levels[i]), n = quarter_lag, default = NA)) %>%
        mutate(!!varname_inter := get(varname_lag) * delta_log_payoff_growth_future) %>%
        mutate(!!varname_inter2 := get(varname_lag) * delta_log_payoff_growth) %>%
        ungroup()
    }
  }
  
  pb <- txtProgressBar(min=1, max=n_groups, style=3)
  progress <- function(n) setTxtProgressBar(pb, n)
  opts <- list(progress=progress)
  
  results <- foreach(i = 1:n_groups, .combine = "rbind", .export = c('fn_identifying_stocks_errors'), .packages = 'tidyverse', .options.snow=opts) %dopar% {
    group       <- group_list[i,1] %>% as.numeric()
    
    df_group <- df %>% filter(.[[group_type]] == group)
    
    T_periods <- df_group %>% count() %>% pull() %>% as.numeric() # Total observations
    T_rolling <- T_periods - window_T + 1        # Number of iterations
    
    output_i <- data.frame() # Empty data frame to store results
    
    for (j in 1:T_rolling) {
      df_window <- df_group %>% slice(j:(j + window_T - 1)) ## Selecting relevant subsample
      
      year_start  <- year(min(df_window$date_fred))
      month_start <- month(min(df_window$date_fred))
      year_end    <- year(max(df_window$date_fred))
      month_end   <- month(max(df_window$date_fred))
      
      recovered <- fn_identifying_stocks_PS(df_window, group, controls_levels_lagged, controls_logs_lagged, controls_levels_inter, controls_logs_inter, unlearnable)
      
      recovered <- recovered %>%
                    mutate(year_start=year_start, month_start=month_start,
                           year_end=year_end, month_end=month_end)
      output_i <- rbind(output_i, recovered)
    }
    output_i
  }
  
  close(pb)
  
  return(results)
}

fn_rolling_stocks_PS2_par <- function(df_input, ivol, df_io_share, controls_levels, controls_logs, window_T, unlearnable, quarter_lag, group_type="permno"){
  group_list <- df_input %>% group_by(across(all_of(group_type))) %>% count()
  n_groups    <- group_list %>% nrow() %>% as.numeric()
  
  df <- df_input
  
  controls_logs_lagged <- c()
  controls_levels_lagged <- c()
  controls_logs_inter <- c()
  controls_levels_inter <- c()
  controls_logs_inter2 <- c()
  controls_levels_inter2 <- c()
  
  if (length(controls_logs) != 0 ) {
    controls_logs_lagged <- vector(length = length(controls_logs))
    for (i in 1:length(controls_logs)) {
      varname_lag <- paste("delta_log_", controls_logs[i], sep = "")
      varname_inter <- paste("inter_log_", controls_logs[i], sep = "")
      varname_inter2 <- paste("inter2_log_", controls_logs[i], sep = "")
      controls_logs_lagged[i] = varname_lag
      controls_logs_inter[i] = varname_inter
      controls_logs_inter2[i] = varname_inter2
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname_lag := log(get(controls_logs[i])) - log(dplyr::lag(get(controls_logs[i]), n = quarter_lag, default = NA))) %>%
        mutate(!!varname_inter := get(varname_lag) * delta_log_payoff_growth_future) %>%
        mutate(!!varname_inter2 := get(varname_lag) * delta_log_payoff_growth) %>%
        ungroup()
    }
  }
  
  if (length(controls_levels) != 0 ) {
    controls_levels_lagged <- vector(length = length(controls_levels))
    for (i in 1:length(controls_levels)) {
      varname_lag <- paste("delta_", controls_levels[i], sep = "")
      varname_inter <- paste("inter_", controls_levels[i], sep = "")
      varname_inter2 <- paste("inter2_", controls_levels[i], sep = "")
      controls_levels_lagged[i] = varname_lag
      controls_levels_inter[i] = varname_inter
      controls_levels_inter2[i] = varname_inter2
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname_lag := get(controls_levels[i]) - dplyr::lag(get(controls_levels[i]), n = quarter_lag, default = NA)) %>%
        mutate(!!varname_inter := get(varname_lag) * delta_log_payoff_growth_future) %>%
        mutate(!!varname_inter2 := get(varname_lag) * delta_log_payoff_growth) %>%
        ungroup()
    }
  }
  
  pb <- txtProgressBar(min=1, max=n_groups, style=3)
  progress <- function(n) setTxtProgressBar(pb, n)
  opts <- list(progress=progress)
  
  results <- foreach(i = 1:n_groups, .combine = "rbind", .export = c('fn_identifying_stocks_errors2'), .packages = 'tidyverse', .options.snow=opts) %dopar% {
    group       <- group_list[i,1] %>% as.numeric()
    
    df_group <- df %>% filter(.[[group_type]] == group)
    
    T_periods <- df_group %>% count() %>% pull() %>% as.numeric() # Total observations
    T_rolling <- T_periods - window_T + 1        # Number of iterations
    
    output_i <- data.frame() # Empty data frame to store results
    
    for (j in 1:T_rolling) {
      df_window <- df_group %>% slice(j:(j + window_T - 1)) ## Selecting relevant subsample
      
      year_start  <- year(min(df_window$date_fred))
      month_start <- month(min(df_window$date_fred))
      year_end    <- year(max(df_window$date_fred))
      month_end   <- month(max(df_window$date_fred))
      
      recovered <- fn_identifying_stocks_errors2(df_window, group, controls_levels_lagged, controls_logs_lagged, controls_levels_inter, controls_logs_inter, controls_levels_inter2, controls_logs_inter2, unlearnable)
      
      recovered <- recovered %>%
        mutate(year_start=year_start, month_start=month_start,
               year_end=year_end, month_end=month_end)
      output_i <- rbind(output_i, recovered)
    }
    output_i
  }
  
  close(pb)
  
  return(results)
}

fn_rolling_stocks_PS_full_par <- function(df_input, ivol, df_io_share, controls_levels, controls_logs, window_T, unlearnable, quarter_lag, group_type="permno"){
  
  group_list <- df_input %>% group_by(across(all_of(group_type))) %>% count()
  n_groups    <- group_list %>% nrow() %>% as.numeric()
  
  n_outputs <- 17
  
  df <- df_input
  
  controls_logs_lagged <- c()
  controls_levels_lagged <- c()
  controls_logs_inter <- c()
  controls_levels_inter <- c()
  
  if (length(controls_logs) != 0 ) {
    controls_logs_lagged <- vector(length = length(controls_logs))
    for (i in 1:length(controls_logs)) {
      varname_lag <- paste("delta_log_", controls_logs[i], sep = "")
      varname_inter <- paste("inter_log_", controls_logs[i], sep = "")
      controls_logs_lagged[i] = varname_lag
      controls_logs_inter[i] = varname_inter
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname_lag := log(get(controls_logs[i])) - log(dplyr::lag(get(controls_logs[i]), n = quarter_lag, default = NA))) %>%
        mutate(!!varname_inter := get(varname_lag) * delta_log_payoff_growth_future) %>%
        ungroup()
    }
  }
  
  if (length(controls_levels) != 0 ) {
    controls_levels_lagged <- vector(length = length(controls_levels))
    for (i in 1:length(controls_levels)) {
      varname_lag <- paste("delta_", controls_levels[i], sep = "")
      varname_inter <- paste("inter_", controls_levels[i], sep = "")
      controls_levels_lagged[i] = varname_lag
      controls_levels_inter[i] = varname_inter
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname_lag := get(controls_levels[i]) - dplyr::lag(get(controls_levels[i]), n = quarter_lag, default = NA)) %>%
        mutate(!!varname_inter := get(varname_lag) * delta_log_payoff_growth_future) %>%
        ungroup()
    }
  }
  
  pb <- txtProgressBar(min=1, max=n_groups, style=3)
  progress <- function(n) setTxtProgressBar(pb, n)
  opts <- list(progress=progress)
  
  results <- foreach(i = 1:n_groups, .combine = "rbind", .export = c('fn_identifying_stocks_PS_full'), .packages = 'tidyverse', .options.snow=opts) %dopar% {
    group       <- group_list[i,1] %>% as.numeric()
    
    df_group <- df %>% filter(.[[group_type]] == group)
    
    T_periods <- df_group %>% count() %>% pull() %>% as.numeric() # Total observations
    T_rolling <- T_periods - window_T + 1        # Number of iterations
    
    output_i <- matrix(ncol = n_outputs, nrow = T_rolling) # Empty data frame to store results
    
    for (j in 1:T_rolling) {
      df_window <- df_group %>% slice(j:(j + window_T - 1)) ## Selecting relevant subsample
      
      year_start  <- year(min(df_window$date_fred))
      month_start <- month(min(df_window$date_fred))
      year_end    <- year(max(df_window$date_fred))
      month_end   <- month(max(df_window$date_fred))
      
      recovered_list <- fn_identifying_stocks_PS_full(df_window, group, controls_levels_lagged, controls_logs_lagged, controls_levels_inter, controls_logs_inter, unlearnable)
      
      output_i[j, 1:4] <- c(year_start, month_start, year_end, month_end)
      output_i[j, 5:n_outputs] <- unlist(recovered_list, use.names = FALSE)
      
    }
    output_i
  }
  
  close(pb)
  
  output <- data.frame(results)
  colnames(output) <- c("year_start", "month_start", "year_end", "month_end", 
                        group_type, "n_obs", "n_factors", "R2_long", "R2_short",
                        "beta_0", "beta_1", "zeta_0", "var_eps", "out_long", "out_short", 
                        "tau_pi", "tau_pi_R")
  
  return(output)
}

fn_rolling_stocks_PS2_full_par <- function(df_input, ivol, df_io_share, controls_levels, controls_logs, window_T, unlearnable, quarter_lag, group_type="permno"){
  
  group_list <- df_input %>% group_by(across(all_of(group_type))) %>% count()
  n_groups    <- group_list %>% nrow() %>% as.numeric()
  
  n_outputs <- 17
  
  df <- df_input
  
  controls_logs_lagged <- c()
  controls_levels_lagged <- c()
  controls_logs_inter <- c()
  controls_levels_inter <- c()
  controls_logs_inter2 <- c()
  controls_levels_inter2 <- c()
  
  if (length(controls_logs) != 0 ) {
    controls_logs_lagged <- vector(length = length(controls_logs))
    for (i in 1:length(controls_logs)) {
      varname_lag <- paste("delta_log_", controls_logs[i], sep = "")
      varname_inter <- paste("inter_log_", controls_logs[i], sep = "")
      varname_inter2 <- paste("inter2_log_", controls_logs[i], sep = "")
      controls_logs_lagged[i] = varname_lag
      controls_logs_inter[i] = varname_inter
      controls_logs_inter2[i] = varname_inter2
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname_lag := log(get(controls_logs[i])) - log(dplyr::lag(get(controls_logs[i]), n = quarter_lag, default = NA))) %>%
        mutate(!!varname_inter := get(varname_lag) * delta_log_payoff_growth_future) %>%
        mutate(!!varname_inter2 := get(varname_lag) * delta_log_payoff_growth) %>%
        ungroup()
    }
  }
  
  if (length(controls_levels) != 0 ) {
    controls_levels_lagged <- vector(length = length(controls_levels))
    for (i in 1:length(controls_levels)) {
      varname_lag <- paste("delta_", controls_levels[i], sep = "")
      varname_inter <- paste("inter_", controls_levels[i], sep = "")
      varname_inter2 <- paste("inter2_", controls_levels[i], sep = "")
      controls_levels_lagged[i] = varname_lag
      controls_levels_inter[i] = varname_inter
      controls_levels_inter2[i] = varname_inter2
      df <- df %>%
        group_by(across(all_of(group_type))) %>%
        mutate(!!varname_lag := get(controls_levels[i]) - dplyr::lag(get(controls_levels[i]), n = quarter_lag, default = NA)) %>%
        mutate(!!varname_inter := get(varname_lag) * delta_log_payoff_growth_future) %>%
        mutate(!!varname_inter2 := get(varname_lag) * delta_log_payoff_growth) %>%
        ungroup()
    }
  }
  
  pb <- txtProgressBar(min=1, max=n_groups, style=3)
  progress <- function(n) setTxtProgressBar(pb, n)
  opts <- list(progress=progress)
  
  results <- foreach(i = 1:n_groups, .combine = "rbind", .export = c('fn_identifying_stocks_PS2_full'), .packages = 'tidyverse', .options.snow=opts) %dopar% {
    group       <- group_list[i,1] %>% as.numeric()
    
    df_group <- df %>% filter(.[[group_type]] == group)
    
    T_periods <- df_group %>% count() %>% pull() %>% as.numeric() # Total observations
    T_rolling <- T_periods - window_T + 1        # Number of iterations
    
    output_i <- matrix(ncol = n_outputs, nrow = T_rolling) # Empty data frame to store results
    
    for (j in 1:T_rolling) {
      df_window <- df_group %>% slice(j:(j + window_T - 1)) ## Selecting relevant subsample
      
      year_start  <- year(min(df_window$date_fred))
      month_start <- month(min(df_window$date_fred))
      year_end    <- year(max(df_window$date_fred))
      month_end   <- month(max(df_window$date_fred))
      
      recovered_list <- fn_identifying_stocks_PS2_full(df_window, group, controls_levels_lagged, controls_logs_lagged, controls_levels_inter, controls_logs_inter, controls_levels_inter2, controls_logs_inter2, unlearnable)
      
      output_i[j, 1:4] <- c(year_start, month_start, year_end, month_end)
      output_i[j, 5:n_outputs] <- unlist(recovered_list, use.names = FALSE)
      
    }
    output_i
  }
  
  close(pb)
  
  output <- data.frame(results)
  colnames(output) <- c("year_start", "month_start", "year_end", "month_end", 
                        group_type, "n_obs", "n_factors", "R2_long", "R2_short",
                        "beta_0", "beta_1", "zeta_0", "var_eps", "out_long", "out_short", 
                        "tau_pi", "tau_pi_R")
  
  return(output)
}